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A COMPUTER PROGRAM VERSION OF THE 
BROUWER ORBITAL THEORY WITH 
OPTIONAL MODIFICATIONS 


1. INTRODUCTION 

This report describes a program which computes osculating elements and 
position-velocity vectors in an earth satellite orbit according to the Brouwer 
theory (Brouwer 1959). Options for certain modifications have been added. 

Basic orbital data are written on a tape called orbital tape which may be used 
as input to other programs for further processing. 

The program written in Fortran IV for the IBM 360 computer is an adapta- 
tion of a program constructed and described by Repass and Chaplick (1965)*. It 
is intended to further revise and extend the program. In view of these intentions 
this report is to be considered as an interim report. No concentrated effort has 
been made to remove all imperfections in the documentation and to check out the 
program thoroughly. This will be done with the final version. Spot checks have 
shown that the program is working correctly and is giving correct results. 

2. THE BROUWER THEORY 

The Brouwer theory (Brouwer 1959) is a first order theory giving expres- 
sions for the osculating values of the Keplerian elements of a satellite orbit. 

To some degree of accuracy the earth's potential at a point at distance r 
from the earth's center of mass and at latitude /3 with respect to the equator is 
given by (IAU 1962) 


u = tL 


r 



£ £)" 
n = 1 


(sin /3) 


(1) 


where P n (x) is the Legendre polynomial of degree n. The constant ft = GM E is 
the product of the gravitational constant G and the mass M E of the earth and has 
the dimension length cubed/time squared. The quantity R is the mean equatorial 
radius of the earth. The coefficient J is called the coefficient of the zonal 

n 

harmonic of order n or briefly the zonal harmonic of order n or the n'th har- 
monic. Since the center of mass of earth is taken as origin J =0. The J n are 
small quantities. It is customary to consider J 2 as \ small quantity of the first 
order and all J n with n > 2 as small quantities of at least second order. 


‘The primary purpose of this program is not to study orbits but to use the Brouwer theory in an 
application. 
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There are also longitude dependent terms in U which Brouwer neglects 
together with the zonal harmonics of order n > 5. For the purposes of his theory 
he writes 


_ m ' uk 2 .. • -> mA 30 


U =H + 

r 


(1-3 sin 2 (3) + 


( 


3 . a 5 . 

sin p + — sui 

2 2 



^ (l - 10 sin 2 /3 + sin 4 /3j 


( 2 ) 


^ A so / 15 . a 35 . 3 - 63 . c A 
+ — — sin /3 - — sin J f3 + — sin 5 /3J 

where 

k 2 = 

A30 = " J3R3 (3) 

k « =-|j,R 4 

A S0 = - Js RS 

In addition to the constant f± we shall also use the constant k defined by 

k 2 = fJL (4) 

If a spherical earth were the only attracting body the potential due to the 
earth would be /+/ r, i.e. all J n with n > 2 would be zero. The orbit of the satel- 
lite would be Keplerian, i.e. the motion would take place in a fixed elliptic orbit 
according to the laws of Kepler. Using well known formulae, the rectangular 
coordinates and velocity components in megameters and megameters per hour 
respectively with respect to the equator and equinox of some epoch may be com- 
puted from the Keplerian elements 

a = semi-major axis in megameters 
e = eccentricity 


2 



t 



I = inclination 

0 = right ascension of ascending 
ci) - argument of perigee 
M = mean anomaly 

The elements a, e, I, are constants whereas M is a linear function of 
the time given by 


node 


I with respect to 
[ equator and equinox 


v in degrees 


M = M 0 + Mj t (5) 

where 

Mj = k a' 3 2 (6) 

is the mean motion of the satellite. 

If r is measured in hours and Mj in degrees per hour then r = 360/M l is 
the period of revolution of the satellite in its closed orbit. 

If a non-spherical earth is assumed, i.e. if the J n with n > 2 are no longer 
all assumed to be zero, the same formulae for computing the coordinates and 
velocity components may be used but the quantities a, e, I , fi , a) occurring in 
these formulae are no longer constants and M is no longer a linear function of 
the time. The variable values 


a, e, I, fi, u), M 

occuring in these formulae under the assumption of a non-spherical earth are 
called the osculating values of the Keplerian elements and are solutions of cer- 
tain differential equations. It is the purpose of the Brouwer theory to find 
approximate expressions for the solutions of these differential equations under 
the assumption that all J n with n > 5 are zero. The deviations of the expres- 
sions for the osculating values from the Keplerian values are of order J 2 but 
it must be mentioned that e occurs as a divisor. 

It should be noted that in the following description of the Brouwer theory tho 
unit of length is the megameter, the unit of time the hour and that angles are ex 
pressed in degrees. These are not the units employed by Brouwer in his paper. 
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In the Brouwer theory an epoch T and six constants 


in megameters 




n. 


Cl) , 


in degrees 


M 


o J 


are chosen. Associated with them a solution of the equations of motion for a, e, 
I,fl ,u) , M of the form 


a B = a 0 


+ 8 a 
+ 8 e 
+ 8 I 

= fi o + n io 1 + sn 


e B = e 0 


= l 0 


( 7 ) 


“b = w 0 + w io * + 8c/J 

M b = M 0 + M 10 t + 8 M 

is constructed. The subscript B has been added to indicate that we deal here 
with approximate expressions based on the Brouwer theory proper. 


In these expressions t is time in hours from the epoch. The quantities 


are constants expressed in degrees per hour and are called the secular motions 
of fl B , a> B , M B . The linear functions 

& = n o + n io 1 


a? = cw 0 + a> 10 t (8) 

M"=M 0+ M 10 t 


4 



t 


are called the secular portions of fi B , &j b ,M h respectively. The secular portions 

of a , c , I are the constants defined by 
B b b J 


c 


H 


- 


r 



( 9 ) 


and thus have no linear terms. 

The secular motion Mj 0 of Mb is approximately equal to the mean motion 
the satellite would have if it were moving in a Keplerian orbit close to the actual 
orbit and is of the order of 180° per hour for a close earth satellite. It is com- 
puted to order J No confusion will arise if we call 360/M 10 the period of the 
satellite even though the actual orbit is not closed and the motion is not periodic. 
But after a period which is of the order of 2 hours the satellite will return to ap- 
proximately the same position in space. This is so since Mi 0 deviates only by 
terms at least of the order J 2 from the value it would have in a nearby Keplerian 
orbit and since Q 10 and cu i0 being at least of order J 2 are also small. Since fl 10 
and oj 10 which are computed to order J* are thus much smaller than M J0 it will 
take the secular portions of and o> B to complete an arc of 360° in a time much 
longer than the period of satellite. In general, times of the order of 100 days are 
required. 

The quantities 

& a in megameters 

5 e 

SI. SO, S a;, SM in degrees 

are series of periodic terms and are called the periodic perturbations of the 
respective elements. They are of order J 2 and are computed to that order 
only. 

The Brouwer theory thus represents the osculating elements as the sums 
of secular portions and periodic perturbations. The periodic perturbations are 
not necessarily equal to zero at the epoch (t = 0) so that 


are not necessarily equal to the osculating values of 

a, e, I, fi, a), M 
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at the epoch. A theory could have been constructed where this would have been 
the case but the periodic perturbations would have been much greater and the 
accuracy of the theory would have suffered. 

For this reason 


are called the Brouwer mean elements for epoch T. However, the first three of 
these quantities are independent of the epoch. The latter three change linearly 
with the epoch. In particular 


a 

e 

I 

n 

u> 

M 


o 

o 

0 

0 

0 

0 


Brouwer mean 
Brouwer mean 
Brouwer mean 
Brouwer mean 
Brouwer mean 
B-ouwer mean 


semi-major axis in megameters 
eccentricity 

inclination > 

right ascension of the node ^ 
argument of perigee 
mean anomaly 


in degrees 


There are two types of periodic terms occurring in the periodic terms. The 
long-period terms are terms whose periods are long compared with the period 
of the satellite while the short-period terms have periods of comparable size 
with the period of the satellite or less. In particular terms whose arguments 
are linear combinations with small integer coefficients of the secular portions 
of and w B are long-period terms. 


We thus write 

8 a = 8 a 

S 

8 e = 8, e + 8 e 

L • 

8 I = 8 l I + 8 b I (10) 

Sfi = 8, fi + 8 fi 

L T s 

8 <j) - 8 , co + 8 co 

L T i 

8 M = 8. M + 8 M 

L * 


where the subscripts L and S indicate long- and short-period terms respectively. 
There are no long-period terms in a. 


6 


* 


While long- period terms are explicitly written down as Fourier series whose 
arguments are linear functions of the time, this is not the case for the short- 
period terms where the true anomaly which is not linear in t is used in the 
arguments. 


The long-period portions of the elements are indicated by primes and are 
defined to be the sums of the secular portions and the long-period perturbations. 
Using (8) and (9) and noting that there are no long-period terms in a we find 


e * = + 5 L e = e o + 5 L e 

I'=I" + S L I,I 0 + 6 l I 

n* + 8 L n = n 0 + fi io 1 + S L n 

60 = °0 + W 10 1 + S L W 

M' = M" + S l M = M 0 4 M 10 t + S l M 


( 11 ) 


The Brouwer expressions for the osculating elements are thus equal to the 
sums of the long-period portions and the short-period perturbations. Thus 


a B = a 0 


+ 8 a = a'+8 a 

S 8 

+ S L e + 5 s e = e' + 5 g e 

,s l i.5.i = i- ,u 
n B = n 0 ( n„t ♦ s L n+ = cr + s t n 


0 =0 

C B 




( 12 ) 


6 J , 


= cj 0 + u> l0 t + h L ci) + 8 s co = + 5 § u) 


M B = M o + M 10 t + 8 l M + 8 g M = M' + 5 s M 


Br'.uwer computes the short-period perturbations by making use of the von 
Zeipel method (von Zeipel 19] 6- 1918). This method furnishes the short-period 
perturbations as functions of the long-period portions. However, Brouwer uses 
in these functions the secular portions e”, I" instead of the long-period portions 
e’,1'. Since he intended to find the short-period terms to order J 2 only and 
since nis procedure causes an error only of order J* this is legitimate. We 
shall come back to this point on page 8. 
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It should be noted that Brouwer does not compute the periodic perturba- 
tions directly but computes first the periodic perturbations in the Delaunay 
variables 


L, G. H, 1, g, h 


which are defined by 


L = (a M ) 1/J 
G = L (1 - e 2 ) 1/2 
H = Geos I 

(13) 

1 - M 

g = oj 

h - Cl 


3. OPTIONAL MODIFICATIONS OF THE BROUWER THEORY 

The program input contains an option which allows the computaton of M 10 
according to Kozai rather than according to Brouwer. This is accomplished by 
adding to the value according to Brouwer, which is of the second order, the 
difference SM 10 of the value according to Kozai which is likewise used as far as 
to order J 2 minus the value according to Brouwer. Mathematical details con- 
cerning this are given in Appendix C (pp. 189-200). 

There is a further option which states whether the long-period terms are to 
be included or not and another option whether the short-period perturbations are 
to be included or not. If they are to be included two further options will be avail- 
able. In one of these the short-period terms are computed, as Brouwer does, with 
the secular portions e" and I", instead of with the long-period portions e' and I' as 
required by the theory. The error (see page 7) is of the second order and is of 
no consequence if one limits oneself as Brouwer does to short-period terms of the 
first order. However, it may be desirable to use second-order short-period terms 
in the semi-major axis a. Then it is essential to have the short-period terms in a 
computed with e' and I'. The second option provides for the computation of the 
short-period terms with e' and I' in a and also e , I, fi , M . 

The second order short-period terms in a are important since the> ftave an 
appreciable effect on the secular motion of M as computed from the Brouwer 
mean elements corresponding to a given position- velocity vector using Brouwer's 
or Kozai's formulation for the mean motion of the mean anomaly. 
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The values of these short-period terms computed by the program are based • 
on expressions given by Kozai (Kozai 1962) and may be optionally included in any 
run with the program being described here. Mathematical details will be found 
in Appendix C (pp. 189-200). 


For the purpose of studying discrepancies with other orbit generators it will 
be useful to have the option of adding certain terms to the Brouwer expressions 
even though with the added terms the expressions may no longer correspond to a 
gravitational solution within the order specified by Brouwer. 


The terms that may be added are linear, quadratic, and cubic terms in the 
time in the semi-major axis, eccentricity, and inclination. The Brouwer theory 
provides for linear terms in time in the right ascension of the node, the argu- 
ment of perigee, and the mean anomaly. Any additional linear terms and quad- 
ratic and cubic terms may be added. Finally, up to 99 sines and cosines may be 
added to the six elements. The arguments of the sines and cosines are arbitrary 
linear functions of the time. The same arguments must be used for each of the 
six elements. The sum of the sines and cosines for each of the six elements will 
be referred to as the supplementary perturbation of this element. 


The optional additions to the Brouwer theory discussed may be written in 
the form 


8 a 

a 

8 e 

a 

5 I 

a 

s n 

a 

S CO 

a 

S M 


a j t + a 2 (0.01 t) 2 + a 3 (0.01 t) 3 + § 2 a + S jup a 

t + e 2 (0.01 t) 2 + e 3 (0.01 t) 3 + 8 jup e 

I x t +I a (0.01 1) 2 + I 3 (0.01 t) 3 + 8 sup I 

dn x t +n 2 (o.oi t) 2 +n 3 (o.oi t ) 3 + s iup n 

dtoj + co ^ (0.01 t) 2 + co 3 (0.01 t) 3 + 5 gup go 

(8M 10 + dM x ) t +M 2 (0.01 1) 2 + M 3 (O.Olt) 3 + 8, up M 


in a 
in e 
in I 
in fi 
in co 
in M 


(14) 
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i 

- 


The quantities 
a y in megameters per hour 
per hour 

I j , dQ j , do> j , dM j in degrees per hour 


define the optional additional linear terms and 
a 2 in megameters per (100 hours) 2 
e 2 per (100 hours) 2 

2 

I j i fij > , M j in degrees per (100 hours) 

and 

a 3 in megameters per (100 hours) 3 
e 3 per (100 hours) 3 

I 3 ,fl 3 ,u> 3 ,M 3 in degrees per (100 hours) 3 


the optional quadratic and cubic terms. For scaling purposes 100 hours instead 
of 1 hour is the unit of time in these terms. 

The quantity 8M 10 is the optional reduction of the mean motion from Brouwer 
to Kozai. The quantity 8 2 a is the sum of the optional second order short-period 
terms in a and 


8 


sup 


a, 


5 

sup 


e. 


8 

sup 


i.s n, 

’ sup ’ 


5 

sup 


W, 


8 

sup 


M 


are the optional supplementary perturbations. 


4. THE BASIC ORBIT GENERATOR FUNCTION 


The first basic orbit generator function consists in computing complete ex- 
pressions for the osculating values of the elements 


1 with respect to 


a = semi-major axis in megameters 
e = eccentricity 
I = inclination 
n = right ascension of ascending node ^equator and equinox 
m = argument of perigee 

M = mean anomaly j 

at t hours from the epoch by means of the formulae 


in degrees 
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a = a B + S a a 

e = e B + b a a 
I = I B + S . I 

n = fi B + s a fi 

a) = oig + 8 a a) 

M = M n + h M 

B a 


(15) 


Here the first terms of the right hand members are given by (12) and correspond 
to the Brouwer theory proper with the added option that the short-period terms 
may be computed with either e",I" or e',I\ The second terms of the right-hand 
members represent the optional additional terms given by (14). Their inclusion 
may cause the osculating values to be no longer consistent with a gravitational 
solution within the order specified by Brouwer. 


To compute the osculating values of the elements for any t the following 
quantities are required: 

Geophysical constants R = mean equatorial radius of the earth in 

megameters 

GM £ = product of the gravitational constant and 
mass of the earth in km cubed/seconds 
squared 

J 2 , J 3 , J 4 , J $ zonal harmonics of orders 2-5 


Epoch (calendar date and time of day) 

Orbital parameters 

Mean Brouwer elements 

a Q semi-major axis in megameters 
e 0 eccentricity 
I inclination 

0 r 

right ascension of ascending node J 

<■ o 0 argument of perigee 

M t mean anomaly 


■\ 

with respect to 
equator and equinox 


> in degrees 


J 
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Additional Secular Coefficients 

Linear coefficients 

rate a x in megameters per hour 

rate e x per hour 

motion I, and additional motions dfipdoij.dMj in degrees per hour 

Quadratic Coefficients 

a 2 in megameters per (100 hours) 2 
e 2 per (100 hours) 2 

1 2 , Qj , oij , M 2 in degrees per (100 hours) 2 

Cubic Coefficients 

a 3 in megameters per (100 hours) 3 
e 3 per (100 hours) 3 

1 3 , ,cj 3 ,M 3 in degrees per (100 hours) 3 

Further it must be indicated whether the reduction of the mean motion of 
the mean anomaly from Brouwer to Kozai, the long- period perturbations, the 
short-period perturbations, the second order short-period terms in a, and the 
supplementary perturbations are to be included. If supplementary perturbations 
are to be included each argument occurring in them is to be defined by giving 
its value in degrees at an epoch to be specified and its change in degrees per 
hour and the coefficients of the cosine and sine of this argument in each of the 
six elements. The epochs to be specified for the arguments may be different 
for each argument and need not coincide with the epoch of the Brouwer mean 
elements. 


Finally, an indication must be given whether the short-period terms are to 
be computed with e",I" or e',1'. 

If expressions based on the Brouwer theory proper are desired the addi- 
tional secular coefficients are to be put equal to zero, i.e. one has to use 

Bj s e t = lj = dflj = dojj = dMj = 0 
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and choose the options not to include the reduction of the mean motion of the mean 
anomaly from Brouwer to Kozai, to include the long- and short-period perturba- 
tions, to compute the latter with e", I", not to include the second order short- period 
terms in a, and not to include any supplementary perturbations. 

The second basic orbit generator function consists in computing for t the 
rectangular equatorial coordinates 

x , y , z in megameters 

and the equatorial velocity components 

x ,y , z in megameters per hour 

This process includes solving Kepler's equation for the eccentric anomaly. 

For convenience sake the radius vector in megameters and the eccentric 
and true anomalies in degrees are included in the output. 

The results are arranged in form of one or several ephemerides. For each 
ephemeris a starting and ending time and a step size must be specified. For 
each ephemeris point 

x »y » z ,x t y » z 

and the osculating values are written on an orbital tape. The format is described 
in Section 9. Information is also written out on the standard output tape at 
specified frequencies. 


5. ALTERNATIVE OPERATIONS OF THE PROGRAM 


The input necessary for the basic orbit generator functions has been described 
in Section 4. However, two alternatives are available. 

* 

Instead of the Brouwer mean elements the osculating values 


a 

e 

I 

n 

CO 

M 


semi-major axis in megameters 

eccentricity 

inclination 

right ascension of the ascending node 
argument of perigee 
mean anomaly 


} with respect to 
equator and equinox 




> in degrees 
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at the epoch may be used in the input. All other quantities of the input must be 
used as described in Section 4. The program then determines by an iteration 
process those Brouwer mean elements which together with the remaining quantities 
of the input would produce the osculating values of the input. After these Brouwer 
mean elements have been determined the program operates as described in 
Section 4. 

As a second alternative, instead of the Brouwer mean elements, the values 
of the rectangular equatorial coordinates 

x , y , z in megameters 

and of the rectangular equatorial velocity components 
x , y , z in megameters per hour 

at the epoch may be used in the input. All other quantities of the input must be 
used as described in Section 4. The program then determines first the values 
of the osculating elements a, e, I, fi, u>, M at the epoch and then, as in the pre- 
ceding case, by an iterative process, those Brouwer mean elements which would 
produce these osculating values. After these have been found the program operates 
as described in Section 4. 


6. UNITS 

The unit of length used in the input described in Section 4 or the alternatives 
described in Section 5 is the megameter. The unit of time is the hour. Angles 
are given in degrees . Velocity components and coefficients of linear terms are 
given per hour and coefficients of the quadratic and cubic terms in t are given 
per (100 hour) 2 and (100 hours) 3 respectively. 

Other orbit generator programs have used other units. An alternate choice 
of units is possible with the present program. So called canonical units may be 
used instead of the units described above. The canonical unit of length (cul) is the 
mean equatorial radius of the earth and the canonical unit of time (cut) is that 
unit which together with the canonical unit of length causes the constant i± in 
expressions (1) and (2) for the earth's potential to be unity. It is 

R 3/2 (GM ) -1/2 seconds 
if R is measured in cm and GM E in cm 3 sec -2 . 
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Coefficients of the linear terms are then expressed per cut and coefficients 
of quadratic and cubic terms per (cut) 2 and (cut) 3 respectively. 

7. BRIEF DESCRIPTION OF THE PROGRAM 

The main program is designated as P73. It requires the following sub- 
routines (SR) and function programs except for library routines. Subroutines 
are called and function programs indicated by * are merely used by name without 
the command CALL. 


ADLH 

BRWR4 

KOMEAN 

SSWTCH 

ALLOT* 

DJUL* 

MAD* 

SUPPO 

ATANQ* 

DMAD* 

PARA 

SUPP1 

ATANZ* 

ELRV 

POLVAL 

TIMC4 

BBRWR 

HMSRZ* 

RHMSZ 

WRT6 

BRWR1 

JULCAL 

RVELZ 

XKEP* 


BRWR2 

The listings for P73 and the subroutines listed above are given in Appendix A 
(pp. 21-172) . An index appears on p. 22. 

All routines except for 

ADLH, BRWR4 , KOMEAN, POLVAL, SSWTCH, SUPPO, SUPP1, WRT6 

were contained, possibly under a slightly different name, in the original Repass 
program package. 

SR ADLH is a slight modification of a program ADDL described by Agreen 
and Fisher (1968). SR BRWR4 is a modification of BRWR2, and KOMEAN, 
POLVAL, SUPPO, SUPP1, and WRT6 have been designed by the author. SSWTCH 
is a routine simulating the sense switches of the IBM 7094 machines. 

In the listings for P73, BRWR4, and all other routines originating from a 
routine of the Repass program package except in the case of DMAD and MAD, the 
following information or parts of it, often appear: 

Purpose (Brief description of function of program or routine) 

Calling sequence to be used 

Input 

Output 

Reference 

Method 

Restrictions 

Accuracy 

Required subprograms 

Timing 

Analysis 
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This information appears also in the listing for BRWR4. Above information in 
many cases is, however, characterized as not available. The list of required 
subprograms refers only to first level subprograms, i.e. a program called by a 
subprogram of some routine is not listed. 

In case of the routines where a list of required subprograms is not given 
the required subprograms may be ascertained from the compilation. 

A brief description of the present program package follows. 

The main program, P73, reads the input except for the ephemeris specifi- 
cations which are read by TIMC4, and the specifications for the supplementary 
perturbations which are read by SUPPO if required. 

The action of PARA depends on which of the three alternative input options 
discussed in Sections 4 and 5 are used. If Brouwer mean elements are used in 
the input PARA determines the values of the osculating elements at the epoch 
by employing the basic orbit generator functions to be described in some detail 
below. If the values of the osculating elements at the epoch are used in the 
input no processing occurs in PARA. Finally, if the position-velocity vector at 
the epoch is used in the input then the values of the osculating elements at the 
epoch are determined by RVELZ. No further processing occurs. 

Thus in every case, after completion of PARA, the values of the osculating 
elements at the epoch will be available. With the help of ELRV the position- 
velocity vector at the epoch is determined. In case the values of the osculating 
elements at epoch or the position-velocity vector at epoch, were used in the 
input, i.e. if the Brouwer mean elements were not used in the input, they are 
now determined. This will be done by BBRWR using an iterative process. If 
the Brouwer mean elements are used in the input then BBRWR is bypassed. In 
every case a point in the program is reached when the Brouwer mean elements 
are known and the basic orbit generator functions discussed in Section 4 may be 
carried out. 

These functions start out with the Brouwer mean elements. The coefficients 
of the secular and long-period terms are computed in BRWR1. This need be 
done just once regardless how many ephemeris points are to be processed. If 
the reduction of the mean motion of the mean anomaly from Brouwer to Kozai 
is to be included it is found from KOMEAN. 

For the epoch and every ephemeris point to be processed, the expressions 
for the long-period terms are evaluated and the short-period perturbations 
are computed. Then the complete values of the osculating elements are formed 
and the position- velocity vector is obtained. All these tasks are accomplished 
by BRWR2 or BRWR4 according to whether e", I" or e', I' are used in the com- 
putation of the short-period perturbations. If the second order short-period 
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terms in a are to be included a call in BRWR2 and BRWR4 to ADLH is necessary.* 
If supplementary perturbations are to be included then a call in BRWR2 and BRWR4 
to SUPP1 is necessary. After the complete values of the osculating elements 
have been formed in BRWR2 or BRWR4 the position-velocity vector is determined 
with the help of ELRV. One step in this procedure consists in solving Kepler's 
equation which requires the use of XKEP. 

The writing of the orbital tape is done in the main program P73 by WRT6. 

Several routines and function programs have not yet been mentioned in this 
brief description. ALLOT reduces an angle modulo 2v . ATANZ and ATANQ find 
an angle from its sine and cosine or its sine and cosine multiplied by a constant 
respectively. DJUL and JULCAL are used in converting from calendar date to 
Julian day and vice versa. DMAD and MAD compute remainders from divisions. 
HSMRZ and RHMSZ convert from hours, minutes, seconds to radians and from 
radians to days, hours, minutes, seconds respectively. Finally, POLVAL com- 
putes polynomials. 


8. DESCRIPTION OF INPUT DECK 

The cards of the input deck are described in the listings in Appendix A of 
the programs where they are read. For convenience sake the relevant informa- 
tion alone is listed in Appendix B (pp. 173-187). 

Several cases may be run. For each case there are cards numbered 1-9 
with several cards for some numbers, in which case distinguishing letters and 
numbers are used. Some of the cards are omitted under special conditions. 
After the last case a number 1 and a number 2 card must follow. The latter 
must have a 0 in cols. 1-3. 

In each case, one or several ephemerides may be computed. Certain data 
for every ephemeris point are written out on a tape called orbital tape, the 
format of which is described in Section 9. 

The quantity KPR in columns 60-65 of card 9 indicates how often the data 
for ephemeris points are to be written on the standard output tape. 

9. ORBITAL TAPES 

The orbital tape is an EBCDIC tape containing double precision data for 
one or several orbits. Each orbit corresponds to one file on the tape. 


In view of the discussion on page 8 it would not be realistic to include the second order short- 
period terms in a and to compute the short-period terms with e", I*. 
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Each file contains a lead record and orbital data for one or several sets of 
equidistant dates. 

The available orbital data are coordinates, velocities, and osculating elements. 


Format 


Lead record: 


NCASE orbit number 

CASE object + orbit number 

(object number is the integral portion, orbit number 
is the decimal portion) 

NTIME = 0 

NPLAN = 0 

NOTHER = 0 

JA(I), I = 1 8 = 0 

INPUT = 0 


ORBJ = O.D 00 
ORB2 = O.D :»0 


Format IX, 13, F8.3, 15, 214, 813, 14, 2F8.3 


For each ephemeris point there are three record 
Record 1: time in minutes from epoch 


Record 2: 


x 

y 

z 


V in megameters 


y f in megameters per hour 
z 

Format IX, 15, 6D21.14 

time in minutes from epoch 

0J) 00 
0J) 00 
0J) 00 


Format IX, 15, 3D23.16 
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Record 3: 


time in minutes from epoch 


t 


a 

e 

I 

n 

U) 




osculating 

value 

of 


> 


M 


J 


semi-major axis in megameters 
eccentricity 

inclination 1 with respect to 

right ascension of ascending node J equator and equinox 
argument of perigee 
mean anomaly 


1 


in 

degrees 


Format IX, 15, 6D21.14 


There are other programs which produce the same type of tape where the 
quantities which are zero in the lead record are not necessarily zero and where 
record 2 of the orbital data contains accelerations instead of zeros. A future 
version of this program will produce these accelerations on the orbital tape. 


Orbital tapes may be input to other programs. One such programs forms 
the differences of orbital data from two orbits which are represented either by 
two files of one orbital tape or two files from two orbital tapes. 


10. DESCRIPTION OF OUTPUT ON STANDARD OUTPUT 

After printing the title page, the input option' from card 2 of input deck are 
listed. A statement concerning units follows. With this information, the value of 
a constant in another orbit generator program, vis SATPOS (Walther and Wales, 
1967) is included. After this information, a list of the values used for the required 
geophysical constants is given. If supplementary perturbations are to be included 
the arguments and coefficients occurring in these aie listed. Then the input orbi- 
tal parameters in the umts used in the input are listed and repeated in the alter- 
nate units. Information concerning the tolerance and maximum number of itera- 
tions allowed in solving Kepler's equation follows. 

A table showing the coefficients of the secular and long-period terms constitutes 
the next group of data. The basis of the calculation of the mean motion of the 
mean anomaly is given. The values of the Brouwer mean elements, osculating 
elements, and the components of the position-velocity vector at the epoch are then 
shown in several units. 

After this ephemeris data are printed. Preceding the data for the individual 
data points information on the limits of the ephemeris and frequency of printing 
is given. For each ephemeris point, in general, sixteen lines of numerical data 
are given in the form of two blocks. The first contains twelve, the second four 
lines. 
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The first line of the first block contains the date in several forms and the 
time from the epoch in minutes. Lines 2-12 give information on various portions 
of the osculating elements. Line 2 gives the gravitational secular portions, i.e. 
the values of the secular portions as derived from the Brouwer theory proper 
with the possible effect of the optional reduction of the mean motion of the mean 
anomaly from Brouwer to Kozai included. The optional secular terms are given 
in lines 3-5 leading to the total secular portions on line 6. Line 7 gives the long- 
period perturbations leading to the long-period portions in line 8. The short-period 
perturbations, second order short- period perturbations, and supplementary 
perturbations follow in lines 9-11. The complete values of the osculating elements 
are then given in line 12, the last line of the first block. 

The first line of the second block lists the radius vector, eccentric and true 
anomalies. The next two lines show the number of iterations required to solve 
Kepler's equation and the last relative correction obtained for the approximation 
to the solution. The final line of the block contains the values of the components 
of the position-velocity vector. 

There are some options to be specified in the input which allow additional 
output to be printed, i.e. detailed output of some phases of the computations. 

These additional outputs are printed if certain fields of columns on card 2 of the 
input deck contain quantities not zero. If a zero is in one of these fields the 
corresponding output will not be printed. 

If the quantity in cols. 19-21 of card 2 of the input deck is not zero inter- 
mediate output from BRWR1 will be printed. It is not labelled. For identification 
the program listing of BRWR1 must be consulted. 

If the quantity in cols. 43-45 of card 2 of the input deck is not zero inter- 
mediate output from KOMEAN will be printed. It is not labeued. For identifica- 
tion the program listing of KOMEAN must be consulted. 

If the second order short-period terms in a are desired and if the quantity 
in cols. 31-33 of card 2 of the input deck is not zero details of the computation 
of the second order short-period terms in a are given for all ephemeris points 
for which data are given in the standard output. There is little labelling in this 
detailed output, which appears between lines 1 and 2 of the first block of twelve 
lines referred to above. Identification should be made with the help of the pro- 
gram listings of ADLH, BRWR2, or BRWR4. 

Finally, if supplementary perturbations are desired and if the quantity in cols. 
37-39 of card 2 of the input deck is not zero intermediate quantities in the compu- 
tation of the supplementary perturbations are printed for every ephemeris point 
for which data are printed in the standard output. This additional output is labelled 
and appears also between lines 1 and 2 of the first block of twelve lines referred 
to above. 
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INDEX 


t 



Page 

P73 

23 

ADLK 

54 

ALLOT 

56 

ATANQ 

60 

ATANZ 

64 

BBRWR 

67 

BRWR1 

72 

BRWR2 

81 

BRWR4 

91 

DJUL 

102 

DMAD 

106 

ELRV 

109 

HMSRZ 

114 

JULCAL 

118 

KOMEAN 

122 

MAD 

130 

PARA 

133 

POLVAL 

138 

RHMSZ 

141 

RVELZ 

145 

SSWTCH 

150 

SUPPO 

153 

SUPP1 

158 

TEMC4 

162 

WRT6 

167 

XKEP 

170 
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APPENDIX B 

DESCRIPTION OF INPUT DECK 
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APPENDIX C 

Mathematical Details of Modification of the Brouwer Theory 


The unit of length in this section is the mean equatorial radius of the earth. 
Let the quantities designated in Section 2 by 


be designated here by 


A", E", I" 

They are the constant terms in the expressions in the Brouwer theory for a , e, 
I . Introduce the quantities 


e", P", 0" 

by means of 

e" = (1 - E" 2 ) 1/2 

P" = A" e " 2 (C-l) 


0" = cos I" 

The Keplerian mean motion n Q is then given by 

n Q = M 1/2 A"' 3/2 (C-2) 

where n is the constant appearing in (1). With the adopted unit for lengths 
and any unit for time m is determined. 

The mean motion of the mean anomaly -£ B as used in the Brouwer theory is 
the first time derivative of the secular portion V B of i B and is given by 

di" 

(C-3) 

i jiCCELXNG r aGC LLnNK. N°T FiLMfeD. 
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where 


B 21 = |p"- 2 e" (-1 + 3©" 2 ) 

B = — P"- 4 «" [(- 15 + 16 e" + 25 e" 2 ) 

128 


+ (30 - 96 e" - 90 e" 2 ) 0" 2 ( c “ 4 ) 

+ (105 + 144 e" + 25 e" 2 )©" 4 ] 

B. = - -i^-P"" 4 e" E" 2 [3 - 30 0" 2 + 35 0" 4 ] 

4 128 


Kozai like Brouwer uses as basic dependent variables the Delaunay variables 

L, G, H, 1, g. h 


defined by (13). 

The complete solutions for L, G, H, 1, g, h according to Kozai are of the 
form 

L = L" + 5L , 1 = 1" + 51 

G = G" + SG , g = g" + 8g (C-5) 

H = H" , h = h" + 8h 


Here L" ,G" , H" are constants and 1" , g" , h" are linear functions of the time. 
Whereas 3G , SI , Sg , Sh contain long- and short-period perturbations, there 
are only short-period terms present in SL. 

The secular portions of a , e , i are designated by double primes and equa- 
tions for the quantities 


are obtained from(13) by replacing L, G , H by L" , G" , H" and a , e , I by a" , 
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Then p" , tj" , 0" are introduced by means of 

p" = a" (1 - e" 2 ) 

V = — = (1 - e" 2 l 1/2 (C-6) 

L" 

9" = cos i" = H"/G" 

The quantity n" is defined by 

n -^ 2 L .- 3 =fl ,/ 2 a - 3/ 2 (C _ ?) 

For the computation of the mean motion of the mean anomaly from equa- 
tion (7.1) in Kozai's article seventeen functions M , i = 0, 1, ... 16 are intro- 
duced. They are defined by 

M n = 1 - 5 9 " 2 Mj = e" (1 - 9 "•*) (1-5 9 u2 )~ 2 = e " (1 - 6>" 2 )'M 2 

M 2 = 3 C 1 - & 2 ) M 3 = 1 - 15 er 2 M 4 = 5 (1 - 7 9- 2 ) 

M s = 2 e" 2 (1 - 5 9" 2 )- 1 = 2 e" 2 /M 0 

M 6 = (l-43d?'' 2 + 155d?" 4 -225d?" 6 )(l - 15d?* 2 ) = (1 _ 43 9" 2 + 155 9“ 4 - 225 6“ 6 ) M 3 
M ? = 10 (1 - 42 9" 2 + 408 9" 4 - 1270 9" 6 + 1575 9“ 8 ) (08) 

M g = 25 (1 - 19 d?" 2 + 75 0" 4 - 105 6 ) (1 - 7 d?" 2 ) = 5 (1 - 19 (?" 2 + 75 4 - 105 <9" 6 ) M 4 

M 9 = e"- 1 (1 - Q" 2 - 2 e" 2 9" 2 ) M 10 = 16 [(32 - 59 9" 2 ) + 6 e" 2 (8 - 13 9" 2 )] 

M u = 4 (1 -9" 2 ) (1 - 5 9" 2 )~ 2 = jM 2 /M 2 M 12 = e" 2 (1 - 5 d?" 2 )" 3 = e " 2 /M 3 

M 13 = (4 - 127 9" 2 + 480 9" 4 - 525 9" 6 ) M 3 

M 14 = 10 (4 - 139 9“ 2 + 1273 9“ 4 - 3805 9" 6 + 3675 9'' 8 ) 

M 15 = 5 (4 - 63 9" 2 + 232 9" 4 - 245 9" 6 ) M 4 

M 16 = - (1 - Q " 2 )- l/2 9 • [2 (l -d?" 2 ) + e" 2 (2 - 3 9 U 2 )] 
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Then the auxiliary constants 


a o* e o- 


1 n’ Pn’ 


defined by equations (6.27) - (6.30) in Kozai's article and the constants 
are obtained by successive approximations from 


1+ — d-ej) 172 0-3 05) J 2 
4 Po 


e 0 = e" + (J 2/1024 a" p" 3 ) M, 



+ M 


5 



(C-9) 


+ (J 2 /] 2 2 ) (1/16 a" p">> M 9 


i 0 = i" - (Jj/ 2048 p" 4 ) sin 2 i" 



+ (J|/J 2 2 ) (1/16 P" 2 ) M 16 


Po = a o O - e o) 
% = (1 - e 2 ) 1/2 


OS 1 , 
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The mean motion of the mean anomaly as used in the Kozai theory will be 
denoted now by d-& K /dt and is equal to the first derivative of the secular portion 
I'l of i and is given according to equation (7.1) in Kozai' s article by 


di" 

7T =n " O + K 21 J 2 + K 22 J| 


K 4 V 


(C-10) 


where 


K 2, -3 tig) 

4 Po 


^22 = — “ — — [10(1-6^5 + 13^ 

128 pj 

-5 eg (5 - 18 + 5 (C-ll) 

- 32 t, 0 (1 - 3 02)2) 

* 4 = " 1287 4 ’ e ° V ° (3 - 30 0 o + 35 *}> 

In the expression for K 4 Kozai's p* has been replaced bvp 4 . 

If the Brouwer theory is modified so that the second order short-period 
terms in a are included, then it may not be unreasonable to assume that 

a" = A" (C-12) 

and that 

e" = E", i" = I" to order J 2 (C-13) 

Then from (C-2), (C-7), (C-12) 

n 0 = n“ (C-14) 
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so that (C-3) and (C-10) lead to 


d d-f 

dt d 


= n 0 {(K 21 - B 21 ) J 2 + (K 22 - B 22 ) J 2 + (K, - B 4 ) J 4 } 


(C— 15) 


According to (C- 13) e" - E" and i" - I" are of order J 2 and according 
to (C-9) the quantities he and hi defined by 


&e=e 0 -e" 
§i = i Q - i" 


(C-16) 


and 


h& = - sin I" S i 


(C-17) 


are of order J 2 . 

It may be conceivable that the differences e" - E" and i" - I" are even 
smaller than the second order of J 2 . Only if we make this assumption will we 
be able to derive a fairly simple approximation for dV K /dt - dV B /dt. The 
assumption implies 


* 0 = E" + he 

i 0 = I* + Si (018) 

d Q = 0" + S0 

where Se and h i are computed from (016) and hd from (017). Using two successive 
iterations of (C9) we find that to order 

a„ = A" [l + yJ 2 -2 7 J jj] (019) 

where 

y = 3 - * (1 - 3 8"*) (020) 

4 P" 2 
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Then we obtain the following correction 


5n = n„ J 


J l (y 3 - 4y B 22 ) - 4 B 4 7 J 2 J 4 



60" 

1 - 30" 2 



(C-21) 


to the mean motion of the mean anomaly computed according to Brouwer. This 
correction will probably make the mean motion to be closer to the value according 
to Kozai. 


To free oneself from above assumptions, one would have to make use of more 
portions of the Kozai theory. This probably would mean so much involvement with 
the Kozai theory that it might just as well be used exclusively and the purpose 
of deriving a reduction from Brouwer's to Kozai's mean motion would be defeated. 

The loi.g period portions L' , G' , H' , V , g' , h' of L , G , H , l , g , h are 
given by 

l' = l" , v = r + b L l 

G' = G" + & l G . g' = g" + & L g (C-22) 

H' -- H" , h' = h' + 8 L h 

there being no long-period perturbations in L and H . 

The quantities 


s l g, s L -t, s L g. S L h 

are the long-period perturbations in G , i , g , h . 
We define 


by means of (13) by replacing L , G , H , a , e , I by the primed quantities L' , G' , 
H* , a' , e' , i' . We also introduce 17 ' , 9‘ by means of 


V = (1 - e' 2 ) 1/2 = G'/L' 
e ' = H"/G* = cos i' 


(C-23) 
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Let f * be the true anomaly and p' = r' /a' the ratio r/a computed from o’ and 
V . Finally define r* , B^ q , B' 2 by 

3u 2 J, 

R* =R* id - 5 f?' 2 ) (f ' -V) 

4 G' 4 

B io = -y(» - 3 < J "> (C-24) 


The perturbations 8L of L , which consist of short- period terms, are then 
given by (see equation (3.8) in Kozai's article) 





b: 


22 



cos 2 (f* 



+ b 2 L (C-25) 


where & 2 L are second order terms due to , J 3 , J 4 . They have been programmed 
for an electronic computer byAgrecnand Fisher (1968). 

The first term of the right hand member of (C-25) represents the first order 
term in &L according to Kozai. Brouwer who uses only first order terms in 8L 
which he then converts to a perturbation in a uses a similar expression but 
instead of g* he uses g' . In order to facilitate comparison of Kozai's expression 
for o L with Brouwer's expression to be given below we rewrite (C-25) and make 
use of (C-24). We then obtain 

8L = SjL + 8'L + 8 2 L (C-26) 


where 

8 L = .— 2 - {(-1 + 3 d' 2 ) (p'- 3 - V" 3 ) + 3(1 - d' 2 ) p'“ 3 cos 2(f' + g')} (C-27) 
4L'a' 

and 


p** J ~ 

b 2 L = 1 — -r — V 3 (1 - 5 e ' 2 ) (1 - d' 2 ) (f' -V) sin 2 (f ' + g') (C-28) 

° r.' 7 _ * 3 
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The quantity S' L represents the reduction of the first order term of (C-25) with 
g* replaced by g' to the form as it appears in (C-25). 

From (13) and (C-26) we obtain to the second order 


a 


1 d t 

«" + Sja + _ 2 L" (o'L + 8 2 L) + _ 


4 a" 


(C-29) 


where 


8 x a = 


2L { (- 1 + 36' 2 ) ( p'~ 3 - V 3 ) + 3(1 - d' 2 ) p'- 3 cos 2 (f * +g')} 


d\ V. f', g*j 

* ( a • 8 ■ ’■ e ) 1 3 82) (t-”' 3 ) 


(C-30) 


+ 3 (1 - d 2 ) !_cos 2 (f + g) 


} 


(C-31) 


Here r /a and f are the ratio of radius vector and semi-major axis and the 
true anomaly computed from the eccentricity and mean anomaly according to 
the laws of Keplerian motion. 

According to the Brouwer theory 


a = A" + S A 

S 


(C-32) 


where 8 s A are the short-period terms in a assumed to be computed with e' , I' . 
This expression is 


5 A = J, 

S J 4 


4> 




0'. 



(C— 33) 



where <t is the same function as the one appearing in (C- 30) , where g' B is the long- 
period portion of g, according to Brouwer, and where r'/A' and f' are the ratio 
of radius vector and semi-major axis and true anomaly computed from the 
long-period portions of the eccentricity and mean anomaly. 

From (C-29) and (C-32) follows 


a" = A" + {(8 s A) 0 - (5j a) Q } 


-I 2 L- [(s;l )o ♦(SjdJ 


4a* 


(C-34) 


where parentheses with subscripts 0 indicate values at the epoch. 

Let us assume that the long-period portions of e , i ,Q , a>, M as computed 
by Brouwer and Kozai differ at most by terms of order J 2> which is a reasonable 
assumption. 

Since in (C-30) and (C-33) occur the same function the difference Sja - S s A is 
at least of order J 2 . Then (C-34) implies that a" - A" will be at least of order 
J 2 . Also, the first two terms in the right hand member of (C-29) will be of the 
same form as the right hand member of (C-32). 

2 

Thus while we cannot prove it, it may be conceivably true that to order J 2 
(C-29) may be replaced by 

1 S 2 A 

a = A" .+ 8 A +-2L" (SiL + 8,L) + -f (C-35) 

M 4 A" 

where 8' L , S 2 L are based on quantities computed according to Brouwer instead 
of Xozai. The terms in (C-35) additional to A" and 8 g A consisting of second order 
short-period terms may lead to an improvement of the Brouwer expression for a. 

We might note again, that the term with S' L allows for the fact that Brouwer 
uses g' while Kozai uses g* in the first order perturbations of L. The term 8 2 L 
represents the second order terms in L in the Kozai theory due to J 2 , J 3 , J 4 
but based on quantities obtained in the Brouwer theory. The last term in (C-35) is 
a second order term due to the conversion from L to a. 

The modifications of the Brouwer theory given by (C-21) and (C-35) will be tested 
experimentally. Results will be reported in a later report. 


200 


A C KNOW LE DG ME NTS 


I wish to acknowledge help received from G. D. Repass, D. Fisher, W. F. 
Cahill, I. J. Cole, and Elise R. Fisher. 


REFERENCES 

Agreen, Russell W., Fisher, David, A Subprogram of the Second-Order Short- 
Period Terms in the Semi- Major Axis, July 1968, GSFC X-552-68-289 

Brouwer, Dirk, Solution of the Problem of Artificial Satellite Theory Without 
Drag, AJ (Astronomical Journal) 64.378.1959. 

IAU (International Astronomical Union), Transactions XI B. 173. 1962 

Kozai, Yoshihide, Second- Order Solution of Artificial Satellite Theory Without 
Air Drag, AJ (Astronomical Journal) 67.446.1962 

Repass, C. D., Chaplick, R. G., A Mutual Visibility Computer Program for 
Communication Satellites, May 1965, GSFC X-547-65-222 

Walter, H. G., Wales, I. M.,The Differential Correction of Close- Earth Satel- 
lite Orbits, Part I. ESRO-SR-7, October 1967. 

v. Zeipel, H., Recherches sur le mouvement des petites planetes, Arkiv for 
Matematik, Astronomi och Fysik 11. Nos. 1, 7; 12. No. 9; 13. No. 3, 
1916-1918. 


NASA - GSFC 


201 



